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Abstract 


It is commonly believed that our own Milky Way is on a collision course with 
the neighbouring Andromeda galaxy. As a result of their merger, predicted in 
around five billion years, the two large spiral galaxies that define the present 
Local Group would form a new elliptical galaxy. Here we consider the latest and 
most accurate observations by the Gaia and Hubble space telescopes, along with 
recent consensus mass estimates to derive possible future scenarios and identify 
the major sources of uncertainty in the evolution of the Local Group over the 
next 10 billion years. We find that the next most massive Local Group member 
galaxies — namely, M33 and the Large Magellanic Cloud — distinctly and radically 
affect the Milky Way - Andromeda orbit. While including M33 increases the 
merger probability, the orbit of the Large Magellanic Cloud runs perpendicular to 
the Milky Way - Andromeda orbit and makes their merger less likely. In the full 
system, we find that uncertainties in the present positions, motions, and masses 
of all galaxies leave room for drastically different outcomes, and a probability of 
close to 50% that there is no Milky Way - Andromeda merger during the next 
10 billion years. 


The Local Group (LG) contains two large spiral galaxies, our own Milky Way (MW) 
and the Andromeda galaxy (hereafter M31), along with approximately 100 known 
smaller galaxies [1]. In addition, it likely hosts additional galaxies yet to be dis- 
covered [2] and, according to the standard cosmological model, a vast number of 
completely dark substructures [3]. 

The negative radial velocity of M31 towards the MW has been known for over a 
century [4], even before its distance was first accurately measured [5]. However, while 
indirect methods had since been used to constrain the transverse components of M31’s 
velocity vector [6, 7, 8], direct measurements of the minute proper motions were only 
achieved much more recently, with the Hubble Space Telescope (HST) [9]. 

The first numerical studies [10] of a possible MW-M31 merger predate even the 
early estimates of the transverse velocity. The finding that the MW-M31 motion 
is close to radial immediately led to the prediction of a likely future collision and 
merger |11, 12, 13]. This scenario has since become the prevalent narrative [14, 15] 
and textbook knowledge [16, 17]. 


Predicting the future of the Local Group 


The MW and M31 both contain remnants of past mergers and interactions with 
other galaxies [18, 19, 20, 21]. Predicting future mergers requires knowledge about the 
present coordinates, velocities, and masses of the systems partaking in the interac- 
tion. In addition to the gravitational force between galaxies, dynamical friction is the 
dominant process in the lead-up to galactic mergers. It describes a transfer of orbital 
kinetic energy to internal energy of the objects involved, and consequently leads to 
the decay of galactic orbits. 

In this study, we parameterise the density profile of each halo as isotropic Navarro- 
Frenk-White profiles [22], use standard assumptions about their concentrations and 
velocity dispersion profiles [23] and calculate dynamical friction using an analytic 
formalism (see Methods for a detailed description). While non-gravitational processes, 
such as gas drag and star formation leading to increased central densities, etc., also 
shape the final stages of a merger, the phase of the orbit that defines the occurrence 
of mergers is largely determined by gravity, which in turn, is dominated by the dark 
matter component in the standard cosmological model [24]. 

The simplest model for the MW-M31 orbit, as considered by [11], contains only 
the two main galaxies. Due to the planar symmetry, only five parameters are required: 
the two masses, along with the initial separation and the two velocity components. In 
general, an orbit with N > 2 galaxies requires specifying 3 x (N — 1) coordinates and 
3x (N—1) velocity components, along with the N masses. More recently, [12] and [13] 
have considered three-body orbits including M33, the third most massive LG galaxy; 
and [13] also considered the Large Magellanic Cloud (LMC), the fourth most massive 
LG galaxy. They still conclude that a merger is certain. We will consider two-body, 
three-body, and four-body systems to study the future evolution of the LG and reveal 
the distinct effects of the M33 and LMC on the MW-M31 orbit. 


A fiducial model based on the most accurate values available 


In the context of the Local Group, it is important to note that apart from the sky 
positions, all parameters including distance moduli, line-of-sight velocities, proper 
motions, masses, and concentrations carry non-negligible uncertainties. We will use 
Monte Carlo (NC) sampling of all values to investigate how these observational errors 
propagate to uncertainties on the future evolution, and in particular, the probability 
of a merger between the MW and M31. 

Our fiducial LG model contains the four most massive LG members: the MW, M31, 
M33, and the LMC, based on the latest and most accurate available data. The mass 
of the MW has recently been extensively studied using Gaia data, with a consensus 
emerging of a total mass close to 1012?Mọo. We adopt a value of Mzo9 = 140.2x10!7Mo 
(excluding the mass of the LMC, which we treat separately). For all other galaxies, 
there is significantly more uncertainty. For M31, we adopt M200 = 1.30.4 x 10!2Mo. 
For M33, we assume M200 = 3+1 x 10!'Mo and for the LMC we assume Mao9 = 
1.5+0.5 x 1014'Mo. See Methods for a review of mass estimates. 

The line-of-sight velocities are well-known and we adopt the values given by [1]. 
For the distance moduli, u, we choose the most accurate and precise recent values 
in the literature: [25] for the LMC and [26] for M33. For M31, we use recent HST 
Cepheid results [27]. 

For the proper motions, where we use the notation, us for the proper motion in dec- 
lination and už, = Hacosô for the proper motion in right ascension, we use the Hubble 
Space Telescope (HST) proper motions of [28] for the LMC, and the combined HST 
and Gaia DR2 proper motions of [13] for M33. Our fiducial M31 proper motions are 
the most precise values in the literature based on Gaia DR3 astrometry [29]. However, 
we find very similar results using the combined Gaia DR2 and HST proper motions of 
[13], as shown in Extended Data Figure 1. All assumed values are listed in Extended 
Table 1. We assume Galactocentric coordinates of RAac = 266.41°, Decgca = —28.94° 
[30], dac = 8.122 kpc [31], and a velocity of (12.9, 245.6, 7.78) kms~! with respect to 
the Sun [30, 31, 32]. 

To account for the fact that the true probability distributions may not be Gaussian 
and to exclude possible effects caused by unrealistic or even unphysical outliers, we 
truncate all distributions at +2ø. This truncation increases the probability of a MW- 
M31 merger by only ~ 10% (see Extended Data Figure 2). For the fiducial model, 
we use 50000 MC samples while for all other variants, we use 2500 samples, ensuring 
that all statistical errors on the merger rate are below 1%. 


The evolution of the MW-M31-M33-LMC system 


The Monte Carlo initial conditions are integrated numerically using a symplectic 
direct scheme (see Methods for details). Figure 1 shows 100 realisations each of the 
MW-M31 orbit, in the two-body MW-M31, and the four-body MW-M31-M33-LMC 
systems. Figure 2 shows the evolution of the distances between the MW and M31 for 
the same sets of orbits, and additionally, for the MW-M31-M33 and MW-M31-LMC 
three-body systems. 


The MW-MB31 two-body orbit evolves in a plane and leads to a merger in slightly 

less than half of cases, the majority of which occur during the second pericentre. The 
addition of M33 increases the merger probability to ~ 2/3, with a similar median 
merger time. However, the addition of the LMC has the opposite effect: the pure MW- 
M31-LMC system experiences a merger in only slightly more than 1/3 of cases and 
the merger probability of the full M31-MW-M33-LMC system is just over 50%. 
For each system, we also show the single orbit obtained by adopting the most 
likely value of each observable, either in the fiducial model or assuming HST+ Gaia 
DR2 proper motions for M31 [13]. In the case of the MW-M31-M33 system, which 
is the one considered by [12] and [13], we reproduce similar orbits and times of first 
pericentre and merger (the differences are likely due to differences in the other model 
parameters). However, despite using newer and more precise measurements, when we 
perform a Monte Carlo analysis, we find considerable uncertainty in the outcome that 
was not previously reported. In particular, in the full MW-M31-M33-LMC system, a 
merger between the MW and M31 occurs within the next 10 Gyr in only approxi- 
mately half of the cases. This is in stark contrast to all previous results that had only 
considered the most likely values without accounting for the numerous and significant 
uncertainties. 

In Figure 3, we show the probability distributions of the merger time and of the 
minimum distance between the MW and M31, as well as the “survival” rate of the 
MW over time, i.e. the probability that no merger with M31 has occurred. In the 
fiducial model, we consider a merger to occur when the distance between any two 
galaxies is below 20 kpc, but we find that our results are not sensitive to this particular 
choice (see Methods for details). Due to the effect of dynamical friction, we find that 
there are two distinct possibilities for the eventual fate of the MW and M31: orbits 
that come within less than ~ 200 kpc eventually merge, which would likely lead to 
the formation of an intermediate-mass elliptical galaxy [10, 33, 34]. For systems that 
merge, we find a median time of 7.6 Gyr in the fiducial model, or 8.0 Gyr adopting 
a 10 kpc threshold. By contrast, orbits with larger pericentres do not decay due to 
dynamical friction. In this case, the MW and M31 continue to evolve in isolation. 
Based on the best current data, both outcomes are almost equally likely. 


The roles of M33 and the LMC 


The distinct effects of each satellite on the MW-M31 orbit are illustrated in Figures 4 
and 5, where we compare the trajectories of the MW and M31 in two-body systems 
and in systems that also include either M33 or the LMC. Both satellites provide 
some additional acceleration in the radial direction of the MW-M31 orbit. However, 
importantly, both satellites also affect the motion of their respective hosts. Including 
M33 in the calculation decreases the transverse velocity of M31 with respect to the 
MW. By contrast, as already pointed out in [35], at its current orbital phase, the 
recoil due to the LMC results in a lower transverse velocity measured between the 
MW and M31. During its short orbital period of ~ 1.5 Gyr, the LMC will accelerate 
the MW to a higher transverse velocity. In addition, the inclusion of M33 largely 
provides momentum in the original MW-M31 plane, while the inclusion of the LMC 
also provides significant momentum perpendicular to the MW-M31 plane. In our 


analysis, the LMC is certain to merge with the MW, and M33 is highly likely (~ 86%) 
to merge with M31 before any possible MW-M31 merger. The net effect of adding 
M33 to the two-body system is to increase the merger probability, while the net effect 
of adding the LMC is to decrease it. 


Sources of uncertainty 


In Figure 6, we show how the merger probability depends on the different observables. 
In each panel, we show the dependence on two variables in the +20 ranges, with the 
remaining variables Monte Carlo sampled. The effects of the concentration parameters 
are shown in Extended Data Figure 3. We find that all masses, the proper motions 
of M31 and M33, and the distance moduli of M31 and M33, significantly impact the 
probability of a merger. The merger probability is positively correlated with the masses 
of the MW, M31, and M33 and negatively correlated with the mass of the LMC. The 
impact of the satellite masses is more pronounced for lower combined masses of the 
MW and M31. The +20 ranges of the M31 proper motions include values that imply 
a merger probability above 90%, but also values that imply a merger probability close 
to zero. The most likely proper motions (assuming no errors) only lead to a merger 
in ~ 2/3 of cases. 

Future more precise proper motion measurements may significantly change the 
expected outcome, although they could make the merger either more or less likely. 
However, if they fall within +1lo of the current most likely values, even precise M31 
proper motion measurements alone will not suffice to determine the outcome. 

Even the comparatively high precision of the line-of-sight velocity and distance 
moduli for M33 and M31 contribute significant uncertainty, with the probability of 
a merger varying between 40% and over 60% for different combinations in the +20 
ranges around the most likely values. 

Given the considerable measurement errors, it is worth noting that cosmologi- 
cal simulations result in a present-day MW-M31 transverse velocity prior of vu, = 
75+% kms~! [36]. There is thus no reason to assume that the transverse velocity mea- 
sured using Gaia DR3 (most likely value of v = 76 kms~*) is overestimated or to 
expect that more precise measurements will result in a lower value. It is also worth 
noting that a perfectly radial present-day M31-MW motion (už, = 38.03 mas yr™t, 
us = —21.37 mas yr~') is neither compatible with current observations nor does it 
result in the highest merger probability in the four-body system. 


Summary 


Even using the latest and most precise available observational data, the future evolu- 
tion of the Local Group is uncertain. Intriguingly, we find an almost equal probability 
for the widely publicised merger scenario (albeit with a later median time to merger) 
and one where the MW and M31 survive unscathed. We reach this conclusion by 
including the LMC and importantly, for the first time, considering the relevant 
uncertainties in the observables. 


Our results are not sensitive to the necessary choices of gravitational softening (see 
Extended Data Figure 4), merger threshold (Extended Data Figure 5), or dynamical 
friction scheme (Extended Data Figure 6)). 

While we have shown that considerable uncertainty results from the proper motion 
measurements of M31, we also find that a more accurate prediction requires more 
precise measurements of the positions, motions, and masses of all participating galax- 
ies. The dependence of the evolution of the MW-M31 system on the treatment of 
other LG galaxies points to further uncertainties. Cosmological simulations suggest 
that ~ 25% of the bound mass of the LG is outside the two main haloes [37]. The 
next most massive individual Local Group galaxies that could impact the MW-M31 
orbit are M32 (an M31 satellite and possible merger remnant [38]) and the Small 
Magellanic Cloud (SMC), a satellite of the MW. Both are at least a factor of five less 
massive than the LMC, and we find that including the SMC, for which proper motion 
measurements are available, does not significantly change the MW-M31 merger prob- 
ability (see Extended Data Figure 7). However, the unaccounted cumulative effects of 
additional substructures as well as that of the cosmic environment introduce further 
uncertainty, particularly towards the far future. An accurate prediction of the evo- 
lution even from perfectly precise observations may require cosmological constrained 
simulations [39, 40, 41] that can account for these effects. Meanwhile, the assumptions 
and simplifications we have made here are likely conservative regarding our central 
claim, that there is considerable uncertainty about the MW-M31 merger. 

Upcoming Gaia data releases will improve the proper motion constraints and mass 
models are continuously refined. However, it is clear that Galactic eschatology is still in 
its infancy and significant work is required before the eventual fate of the Local Group 
can be predicted with any certainty. As it stands, proclamations of the impending 
demise of our Galaxy appear greatly exaggerated. 
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Fig. 1 Possible future MW-M31 orbits. Coloured lines show probability densities for the positions of 
the MW and M31 in 100 Monte Carlo samples of the fiducial model, integrated over 10 Gyr or until 
the merger. On the top and bottom panels, respectively, trajectories are projected in the orbital plane 
and perpendicular to the orbital plane defined by the initial positions and velocities of the MW and 
M31. White markers denote MW-M31 mergers. In the left column, we show the MW-M31 two-body 
system, while in the right, we show the four-body system that includes the MW, M31, M33, and the 
LMC. 
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Fig. 2 Distance between the MW and M31. On each panel, we show 100 realisations of our fiducial 
MC model and also state the probability for a MW-M31 merger within 10 Gyr. In the top row, we 
show the MW-M31 two-body orbits (top-left) and MW-M31-M33-LMC four-body orbits (top-right) 
shown in Figure 1. In the bottom row, we show the MW-M31-M33 (bottom-left) and MW-M31-LMC 
(bottom-right) three-body orbits. White markers denote MW-M31 mergers; percentages indicate the 
fraction of orbits that merge within 10 Gyr. Only slightly more than half of the four-body orbits lead 
to a merger within 10 Gyr. The inclusion of M33 increases the likelihood of a merger, whereas the 
inclusion of the LMC decreases it. White lines show the individual orbits using the most likely values 
of every variable, either assuming the Gaia DR3 proper motions [29] of the fiducial model (solid) or 
HST + Gaia DR2 proper motions [13] (dashed). 
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Fig. 3 Distributions of the merger time, tm, the minimum distance, min(dm31- mw), and the 
“survival” rate of the MW in the MW-M31-M33-LMC system. Blue and red colours show results using 
a distance threshold of 20 kpc (our default) or 10 kpc for a merger, respectively. Mergers happen on 
average later when the distance threshold is lower, but the fraction of systems that merge within 10 
Gyr is similar. The median time to merger is 7.6 Gyr for systems that merge with a 20 kpc threshold 
and 8.0 Gyr for systems that merge with a 10 kpc threshold. The distributions of minimum distance 
show a clear bimodality, independently of the threshold: about half of the systems reach the merger 
threshold within 10 Gyr while the vast majority of the remaining systems do not approach closer than 
200 kpc. The “survival” rate of the MW drops sharply between ~ 6 —9 Gyr and levels off afterwards. 
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Fig. 4 Effects of M33 and the LMC on the trajectory of the MW. As in Figures 1 and 5, panels 
in the top row are projected in the orbital plane defined by the initial MW and M31 positions and 
velocities, while those in the bottom row are projected perpendicular to this orbital plane. Arrows 
point towards the initial position of M31. On the left, we compare MW trajectories in the two-body 
MW-M31 system to those in the three-body MW-M31-M33 system. On the right, we show the orbit 
of the LMC, and that of the MW before and after the merger with the LMC in the MW-M31-LMC 
system. Compared to the two-body orbit, the inclusion of M33 reduces the transverse velocity of 
the MW relative to M31 and introduces only a small velocity perpendicular to the MW-M31 orbital 
plane. By contrast, the LMC increases the MW-M31 transverse velocity and causes a larger velocity 
perpendicular to the MW-M31 orbital plane. 
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Fig. 5 Effects of M33 and the LMC on the trajectory of M31. As in Figures 1 and 4, panels in the 
top row are projected in the orbital plane defined by the initial MW and M31 positions and velocities, 
while those in the bottom row are projected perpendicular to this orbital plane. Arrows point towards 
the initial position of the MW. On the right, we compare M31 trajectories in the two-body MW-M31 
system to those in the three-body MW-M31-LMC system. On the left, we show the orbit of M33, and 
that of M31 before and after a possible merger with M33 in the MW-M31-M33 system. Compared to 
the 2-body orbit, the inclusion of the LMC increases the transverse velocity of M31 and introduces 
motion perpendicular to the initial orbital plane. On the other hand, the inclusion of M33 reduces 
the transverse velocity of M33 with respect to the MW. 
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Fig. 6 Dependence of the merger probability on observables of the MW-M31-M33-LMC system. In 
each panel, we show the probability of a MW-M31 merger within 10 Gyr as a function of two variables, 
with all other observables sampling the probability distributions of the fiducial model. White lines 
on the colour bars indicate the minimum and maximum merger probabilities for the range of values 
shown on each panel, indicating the sensitivity of the merger probability on the two corresponding 
variables. In the top row, we show the dependence on different masses, the second row shows the 
dependence on proper motions, and the third row shows the dependence on distance moduli and 
line of sight velocities. The axes extend to +2ø of the fiducial model. The merger probability is 
positively correlated with the mass of the MW, M31, and M33, and negatively correlated with the 
mass of the LMC. The merger probability strongly depends on us (M31) and už (M31), but also varies 
significantly with us (M33) and už (M33). The uncertainties in the distance moduli for M31 and M33 
also contribute to the uncertainty of the outcome, while the effect of the line-of-sight velocities is 
small. 
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Methods 


Our results are based on numerically integrated initial conditions, which are in turn 
based on Monte Carlo samples of the observational data. We describe here the gen- 
eration of the Monte Carlo samples, the numerical integration, and our treatment of 
dynamical friction and galaxy mergers. We also demonstrate the robustness of our 
findings to the particular choices made and show that our fiducial model and the 
results are conservative in predicting an uncertain future and relatively low merger 
probability. We also present results when the SMC is included in the analysis, in addi- 
tion to the four main galaxies. To facilitate the reproduction of our results and allow 
future work to easily incorporate new observational data, our analysis code is flexible 
and public. 


Monte Carlo Samples 


Our fiducial model consists of the Milky Way and three additional galaxies, M31, 
M33, and the LMC. We assume that the sky coordinates (RA & Dec) for the centres 
of M31, M33, and the LMC are known. We furthermore assume that the position 
of the Galactic centre with respect to the sun is fixed at (Ra, Dec) = (266.4051, - 
28.936175) [30], d = 8.122 kpc [31], and that the galactocentric velocity of the Sun is 
(12.9, 245.6, 7.78) km~? [32, 31, 30}. 

We create Monte Carlo samples for the remaining 20 variables, the four halo 
masses, Mmw, Mm31, Mu33, Mrimc, the four halo concentration parameters, the 
three distance moduli, p, the three sets of proper motions, us and už, and three line- 
of-sight velocities, Vios. Sampling directly in the space of the observables rather than 
sampling in a space of derived variables such as Cartesian coordinates or velocities 
minimises the effect of possible correlations. 

To allow reproducibility and identification of individual orbits across our figures 
when changing parameters of the model such as the set of galaxies included and 
numerical parameters such as the gravitational softening length and merger threshold, 
we re-initialise the pseudo-random number generator with the same values for every 
new set of Monte Carlo samples. While we generally only show the first 100 orbits on 
each plot, all quoted probabilities are computed from at least 2500 samples each, so 
that statistical errors are less than 1%. 

The most likely values of each variable together with the +lo uncertainty are 
either taken directly from single sources in the literature or in the case of masses and 
concentrations (see below), estimated by us based on multiple sources. In creating our 
Monte Carlo samples, we assume that each variable follows a Gaussian probability 
distribution. However, in our fiducial model, we truncate all distributions at +20, 
corresponding to the central ~ 95% of values. While Gaussian distributions are a 
natural assumption for measurement errors, this is not always explicitly stated, and 
it may not reflect the true probability distribution, especially at some distance from 
the most likely values. Indeed, for some variables, untruncated distributions extend 
to unphysical values with finite probabilities. A truncation at 20 also ensures that all 
variables in the samples are physical. 
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Given that our central claim is uncertainty in the future evolution of the LG, 
truncating the probability distributions of the observables is a conservative assump- 
tion that leads to lower uncertainty about the outcome. However, as shown in 
Extended Data Figure 2, our results are not sensitive to the truncation at +20, and 
the merger probability is only slightly lower when the distributions are not truncated. 
On the other hand, even a truncation at only +10 leaves approximately 1/4 of orbits 
that do not merge within 10 Gyr. 


Numerical Integration 


The orbits are integrated using a symplectic leapfrog algorithm in the centre-of-mass 
frame, with a step size of 1 Myr, approximately the time it takes for a galaxy moving 
at 100 kms! to travel 0.1 kpe — 1/200*" of our merger threshold (see below). Our 
results are not affected by the finite time step. 

To account for the fact that the haloes are extended objects, the gravitational force 
between the haloes is softened with a constant softening length of 20 kpc, similar to 
the scale radius of an NFW halo in the mass range we consider here. We also consider 
different choices for the softening and show in Extended Data Figure 4 results with 
softening lengths of 10, 20, or 30 kpc, respectively. A softening length that is too small 
can lead to unphysical hard scattering events during close encounters while a softening 
length that is too large artificially reduces the gravitational force. In the context of 
the Local Group, both of these effects could reduce the merger probability. However, 
we find no strong dependence of the merger probability on the softening length, with 
our adopted fiducial value of 20 kpc resulting in the highest merger probability. 


Dynamical Friction 


To estimate the effect of dynamical friction, we use a modified Chandrasekhar for- 
mula, similar to that used in [24]. The classic Chandrasekhar formula assumes a 
point mass orbiting in the potential of a much more massive host halo, which is in 
turn composed of much less massive particles. This approach has been expanded to 
account for extended satellites [46], and we use the following expression to calculate 
the acceleration of a satellite due to dynamical friction [16] once inside rao99 of the 
host halo: 


2 
dv _ 4rG Mp A sia ee 1 (1) 
dt v? VT v 
where G is the gravitational constant, M is the mass of the satellite, v is the 
velocity of the satellite relative to the host, p is the density of the host at the position 
of the satellite, X = v/(2c) is the ratio between the orbital speed of the satellite and 
the 1D-velocity dispersion, o, of the host at the location of the satellite, and A is the 
Coulomb factor expressed as r/e, with € a scale length that depends on the density of 
the satellite. To determine e€, we adopt an empirical expression derived from N-body 
simulations [47]: 


17 


2.2 rs — 14 kpe ifr, > 8kpc (2) 
E€ = 
0.45 rs if rs < 8kpc 


where r, is the scale radius of the satellite’s NFW halo. Finally, to approximate the 
velocity dispersion of the host at the location of the satellite, we use the expression 
derived in [23] for NFW haloes. 

Standard dynamical friction schemes assume a clear hierarchy between the (much 
more massive) host and the (much less massive) satellite. According to the assump- 
tions underlying Equation 1, the satellite and host enter the calculation in clearly 
defined and distinct roles, with the dynamical friction force applied only on the 
satellite, while the host remains unaffected. 

However, in the Local Group context where galaxies and haloes of similar mass are 
interacting, this introduces an inconsistency. In particular, in the (relatively likely) 
scenario that M31 has a similar mass to the MW, the roles of the satellite and host 
are unclear, but their assignment changes the result of the calculation. For example, 
if the MW is considered the satellite, it would be accelerated by its interaction with 
M31 while M31 would remain unaffected, and only the motion of the MW with respect 
to the other galaxies would be affected while that of M31 would remain unchanged. 
If M31 is considered the satellite, the roles would be reversed. Accelerating only one 
galaxy also violates momentum conservation. 

To make the calculation more symmetrical, in our dynamical friction calculation 
we distribute the dynamical friction force proportionally between the satellite (s) and 
host (h), conserving the total momentum: 


dv, Mnp 


7 My 3 
a” M, 4 Mp’ (3) 
dvp, Ms 

p= 4 
dt aDF M, 4 Mp’ (a 


where v, and M, are the velocity and mass of the satellite, va and Mp are the mass 
and velocity of the host, and apr is the acceleration computed using Equation 1. In 
the limit that the satellite is much less massive than the host, the standard hierarchical 
scheme is recovered and only the satellite is accelerated, while in the limit that both 
galaxies have equal mass, both receive equal and opposite accelerations. 

A small inconsistency remains in that even when the differences in mass are small, 
we still assign the more massive galaxy as the host and the less massive galaxy as the 
satellite when calculating the magnitude of the dynamical friction, where the velocity 
dispersion of the host, ø, but not that of the satellite, and through the Coulomb 
factor, the scale radius of the satellite, but not that of the host, are considered. In our 
spherical halo models, both the velocity dispersion and Coulomb factor depend only 
on the assumed masses and concentrations, and as we discuss below, the concentration 
of the satellite has a greater impact on the dynamical friction force. For two halos 
of significantly different masses, the assignment of host and satellite is clear. For 
an individual case of two haloes of nearly identical masses but different (randomly 
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assigned) concentrations, the choice of calculating the dynamical friction force by 
treating either halo as the satellite seems arbitrary. However, for a large number of 
samples of nearly equal-mass interactions with randomly assigned concentrations, the 
dynamical friction calculations are not biased. We also assume identical distributions 
of concentration parameters for all galaxies. 

In Extended Data Figure 6, we compare the orbits of the fiducial system using 
no dynamical friction, “hierarchical” dynamical friction (only from the more massive 
host to the less massive satellite), and our default “proportional” dynamical friction. 
It is clear that without dynamical friction a MW-M31 merger is highly unlikely. In 
fact, the finite merger rate without dynamical friction depends strongly on our default 
impact parameter threshold of 20 kpc, with a lower threshold, the merger rate can 
become arbitrarily small. On the other hand, when dynamical friction is included, 
the evolution of each orbit is quite similar in the “hierarchical” and “proportional” 
schemes, and the merger rate is not significantly affected by the exact choice of scheme. 

It is worth noting that our semi-analytical approach to dynamical friction is still 
quite simplistic, and while the average behaviour of N-body simulations has been used 
to calibrate parameters, numerical simulations also show that individual systems with 
non-zero internal angular momenta and substructures can have different merger times 
than predicted by these simple equations [48]. A precise prediction of the MW-M31 
orbit will likely require full N-body simulations. On the other hand, we show that even 
a simple dynamical model that assumes no spin, no triaxiality, and no substructure 
results in considerable uncertainty in the future evolution of the MW-M31 system. 


Mergers 


Below a certain distance, the interactions of the gas and stellar components become 
significant, and our simple approach is no longer appropriate for predicting the remain- 
ing orbital evolution. Our aim here is not to predict the precise time of the merger 
(in fact, we argue that such a prediction is futile based on the current data), so we 
simply assume that any system that passes below a threshold distance will eventually 
merge, and identify this time as a lower limit for the time of the merger. 

In our fiducial model, we have adopted a value of 20 kpc as the merger threshold 
for all galaxy interactions. When such a merger occurs between two galaxies, we 
combine the masses and momenta of the two galaxies at their common centre of mass 
and continue the integration. The concentration parameter is set to that of the most 
massive galaxy. 

In Extended Data Figure 5, we compare the sensitivity of our results to adopting 
merger thresholds of 10, 20, and 30 kpc. With a threshold of 10 kpc, we find a slightly 
reduced MW-M31 merger probability. However, raising the threshold from 20 to 30 
kpc has no significant impact on our results. This confirms our assertion in the main 
section: due to the effect of dynamical friction, orbits either inspiral and eventually 
merge, or do not come close enough for dynamical friction to become effective and 
hence do not merge. 

The small fraction of MW-M31 orbits that merge with a threshold of 20 kpc, but 
do not when the threshold is set to 10 kpc, approach with a small impact parameter 
and high velocity. This reduces the effect of dynamical friction and also allows them 
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to escape to a large apocentre. While this may be a possible scenario for the Local 
Group, our methods are not adequate for studying such close interactions between 
galaxies. To be conservative in our prediction of a low merger probability, by setting 
the merger threshold at 20 kpc we assume that these orbits also merge, and with an 
even larger threshold of 30 kpc, we would still predict a similar merger rate. 

For comparison, the best-studied merger of two galaxies with similar stellar masses 
to the MW and M31 are the Antenna galaxies. Their past evolution is reproduced 
with an orbit with pericentre ~ 10 kpc [49], predicted to lead to coalescence ~ 1.3 Gyr 
later [50]. 


Fate of the LMC and M33 


While our main focus here is on the future evolution of the MW-M31 orbit, we 
naturally also make predictions for the evolution of the LMC and M33. 

In the fiducial model, we find that the LMC is certain to merge with the MW 
before any eventual MW-M31 merger. With a merger threshold of 20 kpc, we find a 
median time of the LMC-MW merger of 1.3 Gyr, while with a threshold of 10 kpc, 
we find a median time of 1.9 Gyr. 

For M33, with a merger threshold of 20 kpc, we find an ~ 86% chance of a merger 
with M31 and a median time of 3.3 Gyr, while with a merger threshold of 10 kpc, we 
find an ~ 83% probability of a merger with M31 and a median time of 3.9 Gyr. In 
both cases, we also find a small probability of ~ 1 — 2% for a merger of M33 with the 
MW-M31 remnant after a MW-M31 merger within the next 10 Gyr. 

It must be noted that our simulations are not designed to study these mergers in 
detail and ignore, for example, the impact of the disk of the Milky Way. Nevertheless, 
they broadly agree with the results of [24], who have previously studied the LMC-MW 
encounter in the presence of M31. 


Additional galaxies 


The next most massive LG member galaxy for which proper motion data [28] is 
available is the Small Magellanic Cloud (SMC), a satellite galaxy of the MW that 
has likely been accreted together with the LMC. We repeat our analysis including 
the SMC as a fifth system, and show results in Extended Data Figure 7. Adding the 
SMC, whose mass is approximately 10% of that of the LMC, has no significant effect 
on the merger rate. The SMC properties used are listed along with those of the other 
galaxies in Extended Data Table 1. 


DR2 + HST proper motions 


Due to the fact that the M31 proper motions have the largest impact on the probability 
of a MW-M31 orbit, and for easier comparison with earlier works, particularly [13], 
we also repeat our analysis adopting the HST+Gaia DR2 M31 proper motions of 
[13]. In Extended Data Figure 1, we show the corresponding evolution of the MW- 
M31 distance in the same two-body MW-M31, three-body MW-M31-M33 and MW- 
M31-LMC, and four-body MW-M31-M33-LMC systems. The results can be directly 
compared to Figure 2 which shows the same quantities in our fiducial model that used 
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Gaia DR3 proper motions. In each case, the merger probability is slightly lower using 
the HST+Gaia DR2 proper motions: 40% instead of 44% for the MW-M31 system, 
56% instead of 63% for the MW-M31-M33 system, 34% instead of 37% for the MW- 
M31-LMC system, and 48% instead of 54% for the full MW-M31-M33-LMC system. 
However, both sources of proper motions predict similar distributions of outcomes and 
a similar uncertainty about the MW-M31 merger. This difference can be attributed 
to the slightly lower precision of the HST+Gaia DR2 proper motions. 


Sources for the masses 


As discussed in Figure 6, the assumed masses of all four galaxies and their associated 
uncertainties have a strong impact on the likely evolution of the MW-M31 orbit in 
our fiducial model. Here, we review the mass measurements and show the implications 
for some alternative scenarios. 

MW The total mass of the MW has been extensively studied with different 
methods and tracers, and the accurate astrometry of the Gaia space telescope has 
brought a flurry of recent measurements. Estimates for the total MW mass based on 
Gaia DR2 or DR3 satellite dynamics include 1.171}-75 x 10!2Mo [51], 1.512459 x 
10!?Mo [52], 1.237973 x 101?Mo [53], and 1.14? x 10!2Mo [54] (where the latter 
two works also use a simulation based prior). 

Estimates using rotation curves include 1.08+9:20 x 10!7Mo based on Gaia DR2 
[55], 0.897? 5 x 10!2Mo using stars in the galkin catalogue [56], 0.822+0.052 x 1012?Mo 
using classical Cepheids [57], and 1.0816:1? x 101?Mo from the H3 survey and Gaia 
DR3 [58]. Other recent measurements include 1.547943 x 10!2Mo using combined Gaia 
DR2 and HST kinematics of globular clusters [59], 1.16 + 0.24 x 10'*Mo (including 
the mass of the LMC) using the kinematics of halo stars, 1.267)39 x 10!7M. from 
high velocity RR-Lyrae stars [60], and 1.197339 x 101?Mo from a Bayesian estimate 
using dwarf galaxy kinematics from multiple sources [61]. 

[36] contains an overview of recent measurements and [62] includes a comprehensive 
review of earlier results. 

Analysis from several different tracers, and particularly from the latest studies 
using the latest Gaia observations consistently point towards a Milky Way total mass 
that is close to 10'2Mo. Both the simple mean and median of the above measurements 
are 1.16 x 10!?Mo, but it is worth noting that most methods of measuring the total 
mass of the MW include that of the LMC (at a galactocentric distance of ~ 50 kpc, 
well within fyir or r200), which is not always made explicit. When the LMC is excluded, 
the mass of the MW is reduced by ~ 0.15 x 10'*Mo (see our discussion of the LMC 
mass below). We adopt Myw = 1.00.2 x 10!2Mo excluding the LMC in our fiducial 
model, reflecting the consensus of recent observations. 

M31 While M31 does not enjoy the benefit of accurate Gaia proper motions, 
there nevertheless exist a number of studies estimating its mass. Most recently, [63] 
measured a total mass of 1.14+9:31 x 10!?Mo using rotation curves based on LAMOST 
data release 9 and DESI. [64] measured 1.4 + 0.4 x 10!2Mo derived from satellite 
kinematics, [65] derived 1.2}? x 10!2Mo from kinematics of M31 dwarf spheroidals, 
and [66] measured a total mass of 1.0519:13 x 101?Mo from SED fitting together with 
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the rotation curve and the kinematics of outer globular clusters and satellite galaxies. 
[67] measured 1.0x10!2Mo from the HI rotation curve, [68] measured 2.0703 x 10!7Mo 
from kinematics of the Giant southern stream. [69] measured 1.357972 x 10!2Mo 
and [70] found 1.4*}°3 x 10!?Mo both using outer halo globular clusters, while [71] 
found 0.8+0.1 x 101?Mo from high-velocity planetary nebulae, and [72] found 1.39 + 
0.26 x 101? for the total mass combining disk rotation velocities and radial velocities 
of satellite galaxies and globular clusters. Both the simple mean and median of the 
above values are 1.27 x 10!7Mo. Our adopted mass of Mm31 = 1.3 +0.4 x 10!2Mo in 
the fiducial model reflects the broad consensus of M31 mass estimates using different 
methods, but also the considerable remaining uncertainty. 

M33 Mass estimates of M33 are much more sparse. More than twenty years ago, 
(73, 74] measured a dark matter mass of 5 x 10!°Mo, extrapolated out to a virial mass 
of 5 x 10!!Mo from the measured HI rotation curve, but noted that this results in a 
very low baryon fraction. More recently, [75] obtained a similar result using the Ha 
rotation curve, but noted that because the measurements only extend to a few percent 
of the virial radius, there are no strong constraints on the total dark matter halo. On 
the other hand, abundance matching based on the observed stellar mass results in a 
significantly lower total mass of 1.7 + 0.55 x 10''Mo [76, 77]. Citing both the direct 
measurements and abundance matching, [77] adopt a mass range of 0.8—3.2 x 10'1Mo 
in their dynamical models of the M33-M31 interaction, while [78] and [13] assume a 
total mass of 2.5 x 1014Mo. 

We adopt Mm33 = 3 + 1 x 10!'Mo in the fiducial model, marginally compatible 
with both the extrapolated masses from rotation curve measurements and the results 
of abundance matching and in line with previous studies. However, the results from 
the two methods are certainly in tension. 

LMC For the mass of the LMC, recent mass estimates based on its effect on 
Galactic stellar streams include 1.38793) x 10''Mo [79], 1.30 + 0.3 x 10!1Mo [80], 
1.88733, x 10''Mo [81] and 1.297938 x 10'!Mo. [82] obtain good agreement between 
the perturbations of Milky Way halo stars with an LMC mass of 1.5 x 1014Mo. 
Using the abundance of likely satellites of the LMC, [83] obtain a lower limit of 
1.24 x 10''Mo, and using kinematics of satellites associated with the LMC, [84] find 
1.6573-45 x 10!'M.. Most recently, [85] used 30 LMC globular clusters to infer a total 
mass of 1.807934 x 10'!Mo while [35] find an even higher value of 2.5493 x 104M, 
using a timing argument and to be most consistent with the Hubble Flow around the 
Local Group. The simple mean of the above values is 1.6 x 10!!Mọ while the median 
is 1.5 x 10!!Mọo. A comprehensive recent review on the effect of the LMC on the Milky 
Way, including a discussion of mass measurements, is given by [86], who conclude 
that the LMC mass is likely in the range 1 — 2 x 10''Mo, which matches the choice 
of Mpyo = 1.5 + 0.5 x 10''Mo in our fiducial model. 


Halo Concentrations 


As explained above, we assume that for the purpose of integrating their orbits, all 
galaxies are represented by NF'W haloes [22] with the total masses M209 defined above. 
We assume concentration parameters of 10 + 2 for all galaxies, consistent with results 
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of cosmological simulations [87, 88, 89, 90, 91] in this mass range. While cosmological 
simulations show the average concentration parameter to be mass-dependent, the 
halo-to-halo scatter is significantly larger than the change in the mean concentration 
over this narrow mass range [87]. 

The concentration parameter does not affect the orbital calculation, except for the 
dynamical friction, where (for a given mass), it sets the scale radius, rs = r2o0/c, of 
the “satellite” that enters in Equation 2. 

In Extended Data Figure 3, we show the dependence of the merger probability on 
the concentration. Except for a very low M31 mass, there is only a weak dependence 
of the merger probability on the concentration of M31, which is more likely to be the 
more massive “host” galaxy in the MW-M31 encounter. For the concentration of the 
MW, which, in the MW-M31 interaction is more likely to be the satellite, we find that 
the merger probability is reduced if the concentration is below ~ 8, i.e. below —la 
of our fiducial value. We find no significant dependence on the merger probability on 
the concentration assumed for M33 or the LMC. 
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Extended Data Table 1 Parameters of the Fiducial model 


Galaxy M200 Cc H Ha* Hô Ulos 
[10M] ues rt] [yas ye?) [kma 1] 
MW 100 + 20 102 — = = = 
M31 130 + 40 10+2 24.407 + 0.032 [27] 48.9+10.5 [29] —36.9 + 8.1 [29] —301 +1 [92] 
M33 30 + 10 102 24.67 + 0.07 [26] 31 + 19[13] —29 + 16 [13] —179.2+ 1.8 [1] 
LMC 15+5 10+2 18.477 + 0.026 [25] 1910 + 20 [28] 229 + 47 [28] 262.2 + 3.4 [1] 
SMC 1.5 + 0.5 102 18.99 + 0.03 [93] 722 + 63 [28] —1117 + 61 [28] 145.6 + 0.6 [94] 
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Extended Data Figure 1 Distance between the MW and M31, using HST+Gaia DR2 proper 
motions for M31 [13], analogous to Figure 2. As in Figure 2, solid and dashed white lines denote the 
orbits using the most likely values using either Gaia DR3 proper motions [29] or HST+ Gaia DR2 
proper motions. White markers denote MW-M31 mergers, percentages indicate the fraction of orbits 
that merge within 10 Gyr. The merger rate is slightly lower in all cases when compared to the fiducual 
model that uses the more precise Gata DR3 proper motions, but the results are qualitatively similar. 
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Extended Data Figure 2 Effect of truncating the assumed probability distributions of observables. 
The distance between the MW and M31 (analogous to Figure 2) in the four-body MW-M31-M33- 
LMC system for different truncations of the observables. From left to right, we show results where 
the probability distribution for each observable in the fiducial model is truncated to +lo, +2o (our 
default model, same as in Figure 2), or left untruncated. The fraction of systems that merge is only 
slightly increased when the distributions are clipped at +20 or above. 
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Extended Data Figure 3 Effect of the concentration parameters on the merger probability in the 
four-body MW-M31-M33-LMC system, similar to Figure 6. From left to right, in the top row, we 
show the dependencies on mass and concentration of the MW and M31, while in the bottom row, 
we show those of M33 and the LMC, respectively. The concentration parameter affects the merger 
rate only for the lower mass system in the MW-M31 encounter, which is the MW in most cases. A 
concentration parameter of the MW below ~ 7 (—1.5c), particularly in combination with a low MW 
mass results in a significantly lower merger rate. Unlike their masses, the concentration parameters 
for M33 and the LMC have no discernible effects on the MW-M31 merger probability. 
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Extended Data Figure 4 Effect of the gravitational softening. MW-M31 orbits (analogous to 
Figure 1) and distance between the MW and M31 (analogous to Figure 2) in the four-body MW-M31- 
M33-LMC system for different softening lengths. From left to right, we show results with a softening 
length of 10 kpc, 20 kpc (our default value), and 30 kpc. A softening length that is too small can 
lead to some unrealistically strong kicks in close encounters, while a softening length that is too large 
weakens the overall gravitational attraction. However, the merger fraction is not significantly affected 
by the choice of softening length within this range. 
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Extended Data Figure 5 Effect of the merger threshold on the MW-M31 distance evolution and 
merger rate. The distance between the MW and M31 (analogous to Figure 2) in the four-body MW- 
M31-M33-LMC system for different merger thresholds. From left to right, we show results with a 
threshold of 10 kpc, 20 kpc or 30 kpc for all mergers. The MW-M31 merger probability is not very 
sensitive to the assumed merger threshold, and we obtain almost the same merger probability even 
with a (very generous) threshold of 30 kpc. 
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Extended Data Figure 6 Effect of different schemes of dynamical friction. MW-M31 orbits (topt 
two rows, analogous to Figure 1) and distance between the MW and M31 (bottom row, analogous to 
Figure 2) in the four-body MW-M31-M33-LMC system for different softening lengths. The left column 
assumes no dynamical friction, the middle column uses our default “proportional” scheme where the 
dynamical friction force is divided such that equal and opposite dynamical forces are applied to both 
host and satellite, the right column uses the “hierarchical” scheme where dynamical friction is only 
applied to the less massive object. Dynamical friction is essential for orbits to decay and for the MW- 
M31 merger to occur, but the probability of an MW-M31 merger is not sensitive to the scheme used. 
However, due to its momentum-conserving property, mergers in the proportional scheme are more 
likely to occur close to the original orbital plane compared to the hierarchical scheme. 
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Extended Data Figure 7 The effect of including the SMC on the MW-M31 distance and merger 
rate. In the top row, we show the MW-M31 two-body system, and the MW-M31-M33-LMC four- 
body system in the fiducial model, corresponding to the top row of Figure 2. In the bottom row, 
we add the SMC to both systems, i.e. we show the MW-M31-SMC three-body system and the MW- 
M31-M33-LMC-SMC five-body system. The inclusion of the SMC has only a small effect on most 
MW-M31 orbits and does not significantly affect the total merger rate. 
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